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We study two-species fermion gases with attractive interaction in two-dimensional optical lattices 
producing an array of elongated tube confinements. Focusing on the interplay of Cooper pairing, spin 
imbalance (or magnetization) and intertube tunnehng, we find the pairing gap can exhibit oscillatory 
behavior both along and across the tubes, reminiscent of a Fulde-FerreU-Larkin-Ovchinnikov (FFLO) 
phase. We obtain a Bose-Hubbard-like phase diagram that shows that the magnetization of the 
system undergoes an incompressible-compressible transition as a function of magnetic field and 
intertube tunneling strength. We find the parity of tube-filling imbalance in incompressible states 
is protected by that of the oscillatory pairing gap. Finally, we discuss signatures of this transition 
and thus (indirectly) of the FFLO pairing in cold atom experiments. 
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Experimentally trapping and cooling atoms have led 
to possibilities of studies of condensed quantum many- 
body systems [T]. Tunable interactions via a Feshbach 
resonance [2] have provided realizations of paired super- 
fluid states in ultra cold fermionic gases [3]. Systems 
of attractive majority and minority spin species (f and 
respectively) exhibit various phases distinguishable 
by the Cooper pairing gap, total density, and spin im- 
balance [H [S]. One particularly interesting superfluid 
phase is the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) 
state |6l [7] in which the up spins partially pair with the 
down spins, with a spatially oscillatory pairing ampli- 
tude. However, in three dimensions (3D) the FFLO state 
has remained elusive, while in lower dimensions it can 
hardly be sustained due to the lack of superfluid long- 
range order. Recent attention has turned to an array of 
one-dimensional (ID) systems [5HT^. which is expected 
to have both a wider parameter range for the FFLO (as 
compared with the 3D case) and stabilized long-range 
order through the coupling between these systems. In 
cold atom experiments, such arrays have been realized as 
two-dimensional (2D) optical lattices of elongated tube 
confinements in a global trap [HI [TS] [as illustrated in 
Fig. [l] (a)], with a rich phase diagram observed [151 [H]. 

On the other hand, bosons in optical lattices have 
been extensively studied both theoretically and exper- 
imentally [T]. A well-known model for lattice bosons is 
the Bose-Hubbard (BH) model [T7], which involves tight- 
binding physics in the presence of interactions and yields 
a superfluid-Mott-insulator transition driven by the com- 
petition between localization and itineracy. In this paper, 
we report our findings of a phase diagram of imbalanced 
fermions in lattice tubes that displays a high resemblance 
to the BH phase diagram. Figure [l] (b) shows a contour 
plot of average tube-filling imbalance (number difference 
between the two species in a tube), M, as a function of 
the chemical potential difference (or the magnetic field), 
h, and the intertube tunneling, t. With the decrease in t 
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FIG. 1: (Color online) (a) Illustration of a 2D array of tubes 
in a global 3D trapping potential, filled with f and J, fermions. 
The dashed arrow indicates particle tunneling between tubes, 
while the enclosed pair represents the Cooper pairing, (b) 
Filling imbalance per tube, M, in a radially uniform system 
as a function of chemical potential difference (or magnetic 
field), h, and intertube tunneling, t (obtained from the cor- 
responding Bogoliubov-de Gennes (BdG) calculations). The 
lobe-like plateaux represent magnetic incompressible states 
with corresponding integer fillings, while the outside region 
has finite compressibility, dM/dh. 

across a critical value, the system undergoes a magnetic 
compressible-incompressible (C-IC) transition, with the 
compressibility, dM/dh, going from a finite value to zero, 
accompanied by M changing from fractional to integer 
occupations. In the BH model it is the site-filling num- 
ber that undergoes a transition between C (superfluid) 
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FIG. 2: (Color online) Left: axial profiles of total density p, spin imbalance s (solid light blue and dark red curves, respectively, 
axis on the left of graph), and average magnitude of pairing gap |A| (dashed green curve, axis on the right of graph). Right: 
corresponding profiles of single tubes located in the diagonal from the center (left figure) to the edge (right figure) of the square 
array (represented as in the left, except that the sign of A is revealed). The first and second rows correspond to systems of 25% 
and 50% global polarization, respectively. Top labels show the number, the imbalance and the polarization of the corresponding 
density profiles. All profiles presented are even functions of z. The data obtained were for systems of 2400 particles in a 10 x 10 
tube array with global trapping frequency a; = 0.0625e6//i (which defines the oscillator length, a), tunneling t = 0.014ei,, and 
temperature T = O.left. 



and IC (Mott insulator) states. This similarity implies 
bosonic behavior of the lattice tube fermions, which we 
explain below with a picture of the interplay of the in- 
tratube oscillatory pairing and the intertube coupling. 

The trapped tube lattices filled with spin up or down 
fermions, '0o-=t/^,r(2)j are described by a microscopic 
Hamiltonian 

-t\ ^irV-.r', (1) 

(rr') ,cj 

with the z direction along the tube axis and r = (x, y) 
denoting different tubes (as illustrated in Fig. [l^). The 
one-particle Hamiltonian, = — (/i^/2r7i)9j + m(a;^r^ + 
ujz'^)/2 — /ig., includes the kinetic energy in the z direc- 
tion, the global trapping potential and the chemical po- 
tential for each spin. The on-tube coupling constant is 
given as 5 = 2h^ / [ml'^ {I - 1.033as/f)] in the highly 
elongated tube limit, with being the two-body s-wave 
scattering length and t being the oscillator length of the 
transverse confinement in a tube |18j . The intertube tun- 
neling, given by t = ^Er{Vo/ Enf'^ eM-'^Vyo/ER] 
with Vb being the optical-lattice depth and £'/j being 
the recoil energy (19) . allows particles to move between 
nearest-neighbor tubes, (rr'). 

To analyze the system in the attractive interaction 
regime ((7 < 0), we apply a Bogoliubov-de Gennes (BdG) 
scheme [50] (that has been successfully used in the past to 



describe a variety of tube-confinement problems [9l [2TI - 
I25| and is considered more reliable in the presence of in- 
tertube tunneling), with self-consistent treatment of the 
Hartree potential, Ua — gp-a where pa- is the number 
density, and the BCS pairing field, A = g{ipiil!^). The 
BdG mean-field Hamiltonian is thus written as 

- / E E ^-(^" + U.)^J.r + (A^jrV-lr + H.C.)] 
"'^ r a 

-t [ J2 ^irV-^r'. (2) 

We perform a Bogoliubov transformation, ipariz) = 

J2n{'^riariz)%a - 0-v:^crr{z)li.-a]^ tO rotate Hm tO the 

quasiparticle eigenbasis. We numerically solve the BdG 
equation for self-consistent solutions of the quasiparticle 
wave functions, u and v, as well as their energy spec- 
trum and hence obtain the spatial profiles for total den- 
sity, p = p^ + pi, spin imbalance (or magnetization), 
s — p-f — pi, and pairing gap, A. The following results 
were obtained for an interaction strength set by choosing 
the binding energy £{, = mg^ /Ah? = 16huj. This setup 
can describe a realizable ^Li system of w = 27r x 200 Hz, 
Tc - 100 nK, and 10-100 particles per tube [H]. 

We first look at the spin imbalance distribution and 
the appearance of oscillatory pairing in an isotropic trap 
{ujr — oj). In Fig. [2| we plot the axial profiles of p, s 
and the average of |A| (left side of the figure) of the sys- 
tem after tracing out the r degree of freedom; as well 
as the corresponding profiles in single tubes aligned in 
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diagonal from the center to the edge of the array (right 
side). The first and second rows correspond to a lower 
global polarization (LP) of 25% and a higher one (HP) 
of 50%, respectively. Here the polarization P is the ra- 
tio of the global imbalance, Mg, to the global number, 
Ng. We find that the imbalanced regions of tubes al- 
ways accommodate oscillatory pairing (FFLO), with the 
concurrence of imbalance local maxima and gap nodes, 
except for the gap decaying to zero in entering a fully po- 
larized region {p = s). In the LP case, single-tube profiles 
in the z direction can exhibit an FFLO-pairing center, a 
BCS-like fully paired off-center region (s ~ and no gap 
nodes), a small re-entrance to FFLO, and fully polarized 
tails. A similar trend is seen in the f direction: center 
tubes with a larger FFLO region, off-center tubes with 
a larger BCS-like region, and fully polarized edge tubes. 
The oscillatory pairing behavior is also revealed by the 
sign changes at fixed z across the tubes. In the HP case, 
the structure of an FFLO center and fully polarized tails 
remains, except the BCS-like regions disappear. 

The axial profiles of both the spin imbalance and the 
pairing gap in the HP case exhibit a clearer oscillation 
than those in the LP case. This indicates an align- 
ment of the gap nodes (occupied by unpaired majority 
spins) across the tubes and hence implies a significant 
role of the intertube coupling (via the particle tunnel- 
ing in our model). To investigate the interplay between 
these parameters, we compute the total imbalance per 
tube, M (which also represents the occupation number 
of unpaired majority) of a radially homogeneous system 
{ujr = 0) at zero temperature. Figure [T] (b) shows a di- 
agram of M as a function of intertube tunneling, i, and 
chemical potential difference, h. It exhibits a remarkable 
similarity to the phase diagram of the BH model [TT], 
which describes a lattice system of interacting bosons. 
A salient feature of the BH model is the existence of 
a quantum phase transition between a incompressible 
Mott insulator with integer fillings and a compressible 
Bose-Einstein condensate with factional site occupations, 
driven by an energetic competition between tunneling 
and interaction. Our tube system, in analogy, possesses 
a magnetic C-IC transition, with the two phases identi- 
fied by the filling of imbalance and the compressibility, 
dM/dh. At small tunneling, an increase of the magnetic 
field, h, moves the system from an IC state [one of the 
lobe-like plateaux in Fig. [l] (b)] to another one with dif- 
ferent filling integer, across a narrow C region. If we 
begin with a fixed filling integer, with an increase in tun- 
neling the corresponding plateau shrinks until a critical 
tunneling (the tip of the lobe) , beyond which the system 
enters the C regime; with still the same integer filling. 

Figure [3] shows the M-h curves at several tunneling 
strengths [the corresponding section profiles of Fig. [T] (b)] 
in a wide range of h. We see that in the low-/i regime the 
whole system is a BCS superfluid (M = 0) and the FFLO 
sets in at a critical h beyond which M becomes finite. In 




FIG. 3: (Color online) Imbalance per tube, M, vs magnetic 
field, h at various tunneling strengths. Curves from bottom 
to top correspond to t/et, = 0.003, 0.014, 0.02 and 0.05 (or 
lattice depth Vq/Er = 12, 7, 6 and 3.8, respectively), with 
the top three curves shifted up. Inset: dM/dh vs. h a,t t = 
O.OSeb- A A-like transition appears around P = 17.5%. At 
small t the curves display the same trend but with staggered 
substructures. The data obtained were for systems of 2560 
particles in an 8 x 8 array of tubes at zero temperature. 



the relatively low polarization regime (P = 10%-20%, 
around h = O.Tfiw), the compressibility displays a A-like 
transition with a peak corresponding to P = 17.5% (see 
inset), which agrees with the critical polarization for the 
disappearance of BCS-like outer regions in the experi- 
ment |15j . Our data show that this critical polarization 
is insensitive to the strength of intertube tunneling. 

Let us discuss a difference between our system and the 
BH model: the filling integer between two neighboring 
plateaux jumps by two in our system whereas it jumps 
by one in the BH model. Since the trapping potential 
as an even function of z conserves the parity symmetry 
of the gap profile fM], there are even (odd) nodal solu- 
tions capable of accommodating an even (odd) number 
of unpaired majority spins per tube in an incompressible 
state. Although the results presented in this paper are 
for the even solutions, we find also the odd ones exhibit- 
ing the same BH-like diagram. In either case, the integer 
filling jumps by two. 

Moreover, we find the system prefers all the tubes of 
the same parity, otherwise an energy cost will arise due 
to the formation of domain walls between tubes of op- 
posite parities. Such energy cost not only prohibits the 
system from globally "flipping" between opposite pari- 
ties but also suppresses the single-particle tunneling of 
the unpaired majority spins. Therefore, the leading be- 
havior remaining is a two-body tunneling that conserves 
the tube parity. This bosonic nature allows us to write 
down an effective BH Hamiltonian for the system, 

Hcs = -i ^i^r' + -^(^r) - hM,, (3) 

<rr'> r 
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FIG. 4: (Color online) Response of imbalance (magnetization) 
axial profile to tunneling modulation around ujm ~ ^- Curves 
from bottom to top correspond to t/et = 0.003, 0.014, 0.02 
and 0.05, with each of the top three curves shifted up (as 
illustrated in Fig. |3|. A large (small) response indicates a C 
(IC) state. 

where w creates two unpaired majority spins while pre- 
serving parity and M = Wh is the imbalance operator of 
a fixed parity. The effective tunneling, is estimated as 
the ratio of to the domain wall energy. The on-tube 
energy E can be determined through the BdG calcula- 
tions or be estimated as the trapping energy plus the 
Hartree interaction (the pairing energy is negligible). Us- 
ing Thomas-Fermi density profiles, we obtain the leading 
dependence of E on 

0.025 lnP)]M2 -I- 0{M^). Negative g leads to repulsion 
between unpaired majority spins. Fitting the numerical 
data in Fig. [l] (b) with the mean-field solution of the 
BH model in the small tunneling regime |26j . we find 
E = O.OSfbM^ and i ~ 0.024t For a uniform tube, 
analytic limits of the energy take the form of non lin- 
ear functions in M [57], so we expect the C-IC tran- 
sition in such a tube array too. Because the interac- 
tion plays a key role in the phase transition in the BH 
model, this magnetic C-IC transition is different from 
the commensurate-incommensurate transition identified 
in Ref. [9J that is interpreted as band filling of unpaired 
majority spins. These two different scenarios for the 
gapped states, analogous to the difference between the 
band insulator (cf. Ref. [9]) and the Mott insulator (here), 
could come from differences in the two setups considered. 
For example, a setup of homogeneous tubes and periodic 
boundaries in Ref. [9^ may favor a pairing function in the 
plane wave basis (FF type), while a trap confinement 
in our setup naturally conserves the parity symmetry of 
the pairing in real space (LO type). In addition, the pa- 
rameter regimes considered in these two studies are also 
different. 

We now turn to the issue of probing such transitions 
in experiments. A salient characteristic of an IC state 
is that the energy spectrum exhibits a gap between the 
ground-state energy and the low excitations, which is not 



present in a C state. A realizable technique for explor- 
ing the excitation spectrum is to detect the response of 
the lattice system to a time periodic modulation of the 
lattice depth. This probe has been successfully applied 
to distinguish an IC state from a C one in cold-atom ex- 
periments of quasi- ID lattice bosons [5S] and 3D lattice 
fermions 

For our lattice of tubes, we calculate the linear re- 
sponse of the imbalance in the "hidden" dimension of the 
2D array — the axial profile, 5s{z) — to a modulation in 
tunneling strength, which can be induced by modulating 
the lattice depth. We obtain 5s{z^ujm) in the frequency 
domain as 

n<0,ri'>0 r 

where e„ are quasiparticle eigenenergies and Wm is 
the modulation frequency. The amplitude A^n' oc 

Z]<rr'> ("n'ir'"'4r ~ "n'tr'""Tr)' Combined with the 
orthogonality between quasiparticle wave functions, 
guarantees that 5s{z) has even parity in z and conserves 
the total imbalance (/^ 5s = 0). 

To refiect the natural inhomogeneity in experiments, 
we present our results for a trapped setup, as in Fig. [2] 
of 50% polarization and at a low temperature (T = 
O.Olefc). Figure |4] shows 5s at four different tunneling 
strengths, those illustrated for a radially homogeneous 
case in Fig. |3] It can be seen that in the IC regime 
where the M-h curves display a staggered structure, the 
imbalance profile has little response {5s 0), while in 
the C regime 5s shows larger features signaling a strong 
response. Notice that our results are robust against fluc- 
tuations in total numbers coming from the Josephson 
coupling between tubes, due to the nature of the BdG 
Hamiltonian that conserves the total spin imbalance. 

In conclusion, we identified an interesting transition in 
a lattice array of tubes with a strong analogy to the ex- 
tensively studied compressible-incompressible transition 
in the BH model. In our fermionic system, the underly- 
ing reason for the effective bosonic behavior is crucially 
linked to the oscillatory nature of the pairing amplitude. 
The observation of this transition would thus constitute 
additional evidence of the FFLO state hinted at by recent 
experiments. Our prediction of easily accessible experi- 
mental signatures is quite encouraging and this transition 
deserves further theoretical and experimental investiga- 
tion. 
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